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We study theoretically the influence of the underlying domain microstructure on the electromechanical prop- 
erties of ferroelectrics. Our calculations are based on a continuum approach that incorporates the long-range 
elastic and electrostatic interactions. The theory is used to simulate the piezoelectric properties of a two di- 
mensional model ferroelectric crystal. Simulation results indicate that the electromechanical response of the 
ferroelectric is strongly dependent on the domain microsctructure, including domain walls. This is particularly 
true for the case when an electric field is applied along a non-polar direction. 
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Ferroelectrics are the best piezoelectric materials that can 
convert electrical energy into mechanical energy and vice 
versat This electromechanical property arises due to the cou- 
pling of the spontaneous polarization with lattice strain. Many 
devices such as ultrasonic transducers and piezoelectric actu- 
ators make use of this property 2 . Recently, there has been 
considerable interest in this field due to the observation of a 
giant piezoelectric response when the applied field is along 
a non-polar directional. It is believed that this "superpiezo- 
electric" response is due to the symmetry change caused by a 
rotation of the polarization towards the direction of the applied 
fielc£. Domain configurations produced by the non-polar di- 
rection field are termed engineered domains. There are also 
a large number of domain walls between the degenerate vari- 
ants which affect the piezoelectric property. It is important to 
understand the role played by such domain microstructures on 
the piezoelectric response. 

Electromechanical properties of ferroelectrics have been 
studied theoretically using first-principle calculations^^. A 
continuum Landau theory describing a single domain or ho- 
mogeneous state has been used to study the electromechanical 
properties of BaTiOj, as function of temperature and electric 
field direction^. Although such calculations provide valuable 
insights into the physics of the polarization-strain coupling, 
they do not describe inhomogeneities due to domains and do- 
main walls. In this Letter, we study the behavior of domain 
patterns under applied electric field and investigate how the 
microstructural evolution influences the average electrome- 
chanical response of ferroelectric materials. 

Our approach is to use a time-dependent Ginzburg-Landau 
model2i£ with long-range elastic and electrostatic effects. We 
restrict ourselves to a 2D ferroelectric transition to illustrate 
the basic principles and use parameters from a model for 
BaTiOj, in our calculations^. We find that the electromechan- 
ical response is highly orientation and microstructure depen- 
dent. When the field is applied along one of the polar direc- 
tions, domain switching results in higher strains compared to 
the single domain state. For fields applied along a non-polar 
direction, the domain walls serve as nucleation sites for a field 
induced structural phase transition. 
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where P x and P y are the polarization components. The free 
energy coefficients a\, an, am determine the ferroelec- 
tric phase and the gradient coefficients g±, g 2 and g 3 are a 
measure of domain wall energies. E x and E y are the compo- 
nents of an external electric field. Elastic properties are stud- 
ied by using the strains 771 = r\ xx + g yy , 772 = i] xx - r\ yv 
and 773 = i] xy , where 77^ is the linearized strain tensor de- 
fined as rjij — (ui,j + Uj,i)/2 = x,y), u; being the 
components of the displacement vector. The electromechani- 
cal coupling is described in terms of these strain variables as 
the free energy F em = A / dr[{i h - Qi(P x 2 + P y 2 )} 2 + 

{772 - Q2{P X 2 - Py 2 )} 2 + {% - QsPxPy} 2 }- Here Qi, 
Q 2 and Q3 are obtained from the electrostrictive constants of 
the material as Qi = Q n + Qi 2 , Q 2 = Q n - Q 12 and 
Q3 = Q44 (electrostrictive constants describe coupling be- 
tween strains and polarization as i} xx = Q\\P X 2 + Q\ 2 P y 2 , 
Vyy = QnPy 2 + QnPx 2 and ij xy = Q 44: P x Py). Notice 
that the free energy F em vanishes for a homogeneous state as 
the homogeneous strains in equilibrium are given by 771 e = 

Ql(P X 2 + Py 2 ), V2 C - Q2(P X 2 ~ Py 2 ), and 77 3 G = Q 3 P x Py. 

However, this free energy does not vanish for an inhomoge- 
neous state. For an inhomogeneous state, the strains 771, 772 
and 773 are related to each other by the elastic compatibility 
constraintii V 2 ?7i - ( - 7^)^72 - afaTj^ = 0. Using this 
relation, the strain 771 can be eliminated from F em resulting 
in a nonlocal interaction between the strains involving 772 and 
773. Using the equilibrium strains defined by 772 e and 773 e , the 
electromechanical free energy can be written as 



energy* that describes the ferroelectric transformation and is 
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where the k = mode has been excluded from the above inte- 
gral. The constant A is the strength of this nonlocal interaction 
and hence it influences the underlying microstructure. The 
quantities Ti(k), ^(fc) and Ts(k) are respectively the Fourier 
transforms of Qi{P x 2 + P y 2 ) ,Q 2 (P X 2 -P y 2 ) and Q 3 P x P y , 
C 2 = {k x 2 - k y 2 )/(k x 2 + k y 2 ) and C 3 = k x k y /(k x 2 + k 2 ) 
are the orientation dependent kernels. The electrostatic contri- 
bution to the free energy is calculated by considering the depo- 
larization energyi^ F es =-fJ,J df{E d ■ P + e (E d • E d /2)}, 
where E d is the internal depolarization field due to the dipoles 
and /i is the strength of this interaction. The field E d can be 
calculated from an underlying potential as E d = — V(f>. If 
we assume that there is no free charge in the system, then 
V • D = 0, where D is the electric displacement vector defined 
by D = e^E d + P. This equation gives rise to the constraint 
— eo V 2 + V • P — 0. The potential <fi is eliminated from the 
free energy F es using the above constraint to express F es in 
Fourier space as 

F es = {fx/2e ) J dk\k x P x (k) + k v P y {k)\ 2 . (3) 

The above integral excludes the homogeneous k = mode 
which means that the homogeneous depolarization field due 
to surface charges has been neglected. The total energy is 
defined as F = Fi + F em + F es with two additional constants, 
i.e. A and \i are essential for the description of multi-domain 
states. 

The dynamics of the polarization fields is given by the re- 
laxational time-dependent Ginzburg-Landau equations = 
—Jjp-, where 7 is a dissipation coefficient and i = x,y repre- 
sents the polarization components. We first introduce rescaled 
variables defined as u = P x /Pq, v = P v /Pq, £ = r/5 and 
t* = 7|ai(Tb)|f, where To is a fixed temperature. In this 
work, we use the parameters! for BaTiO^ for the local part 
of the free energy Fi. The parameters which can be depen- 
dent on the temperature T are: a x = 3.34 x 10 5 (T - 381) 
VmCT \ a u = 4.69 x 10 6 (T - 393) - 2.02 x 10 8 Vm 5 CT 3 , 
a ln = -5.52 x 10 7 (T - 393) + 2.76 x 10 9 Vm 9 C" 5 , a 12 = 
3.23 x 10 8 Vm 5 C" 3 and a lu = 4.47 x 10 9 Vm 9 C" 5 . The 
electrostrictive constants are given as Qn = 0.11 m 4 C -2 , 
Q12 = -0.045 m 4 CT 2 and Q 4 4 = 0.029 m 4 C" 2 . We assume 
that the coefficients g\ = g 2 = .93 = g and use the value 
g = 0.025 x 10~ 7 Vm 3 /C quoted in the literature 1 ^. To calcu- 
late the rescaled quantities, we use To = 298-fr' and To = 0.26 
Cm~ 2 and 5 ~ 6.7 nm. The values chosen for the long-range 
parameters are A = 0.25|ai(T )|/T 2 and /i = 20e |ai(T )|. 

The time-dependent Ginzburg-Landau model with the 
above rescaled parameters is used to simulate the domain pat- 
terns and electromechanical properties. The equations are dis- 
cretized on a 128 x 128 grid with the Euler scheme using pe- 
riodic boundary conditions. For the length rescaling factor 
5 ~ 6.7 nm, this discretization corresponds to a system of 
size ~ 0.85 /im x0.85 fim. We first simulate the properties 
of this 2T> model at T = 298A'. At this temperature, the 
minima of the free energy T; define a rectangular ferroelec- 
tric phase with the four degenerate states (±0.26, 0) Cm -2 



and (0, ±0.26) Cm -2 . The dynamical equations are solved 
starting from small amplitude random initial conditions cor- 
responding to a quenched paraelectric phase. Domains of the 
four degenerate states form and a domain growth process takes 
place. Eventually, the growth stops and a stable multi-domain 
state shown in Fig. 1(a) is obtained. All four polarization 
variants exist in this state. Note that there are only 90° walls 
(domain walls across which the polarization angle changes by 
ninety degree) and these walls are aligned along the [11] or 
[11] directions. The domain wall orientations are governed 
by the underlyling symmetry encoded in the anisotropic ker- 
nel in F es . Another interesting feature of this domain struc- 
ture is that there is no charge accumulation. This is clear from 
the fact that there are no head-head or tail-tail configurations 
at the domain walls. At some of the domain junctions polar- 
ization vortices are observed and at other junctions the heads 
coming in appear to balance the tails going out, thereby main- 
taining zero net charge. 

To simulate the effect of an external electric field, the evo- 
lution equations are solved with a varying E y (electric field in 
the [01] direction) while E x is kept at zero. The domain con- 
figuration of Fig. 1(a) is the initial condition as E y is varied 
quasi-statically (the system is allowed to relax for t* = 1000 
after each change) from E y = to E y = 72.07 kV/cm in 
steps of 1.85 kV/cm. It is clear from figures 1(b)- 1(d) that 
this electric field causes the unfavorable domains to switch 
and since the field is along one of the polarization directions, 
eventually (E y ~ 9 kV/cm) a single domain polarized along 
the [01] direction is established. On removing the field, the 
system stays in this single domain state and there is an under- 
lying hysteresis. Figure 2(a) shows the variation of (P x ) and 
(Pa) with E y for the situation shown in Fig. 1 . Also shown is 
the effect of electric field E y on P x and P y for an initial single 
domain with P x = and P y = 0.26 Cm~ 2 (these single do- 
main response curves are obtained by minimizing the homo- 
geneous local free energy T; only). Beyond E y ~ 9 kV/cm, 
the single and multi-domain responses coincide as both corre- 
spond to single domain states. 

The appearance of (P x ) during the transient indicates that 
the polarization reversal of the P y < regions is via rota- 
tion of dipoles instead of direct 180° flipping. Because the 
simulation size is limited, the number of clockwise and coun- 
terclockwise rotating dipoles is different, so that a transient 
(P x ) is observed. This can only be observed in very small 
systems in reality. To study the electromechanical behavior, 
we have also computed the behavior of uniaxial strain with 
the applied electric field. The uniaxial strain is calculated as 

(Vyy( E y)) ~ (%»(0)}» where Vyy = QnPy 2 + QuP x 2 - Fig- 
ure 2(b) shows the variation of the uniaxial strain with the 
electric field. The electromechanical response of the single 
domain state is also shown. It is clear that the multi-domain 
state generates much higher strain compared to the single do- 
main case. The switching of 90° domains provides the ex- 
tra strain in the multi-domain state. The piezoelectric con- 
stant c?33 (calculated as the slope of the strain vs. electric 
field curve) for the large field fully polarized state is nearly 
the same for both the single-domain and multi-domain re- 
sponses (~ 82 pC/N). We have also calculated the dielectric 
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constant e yy = 1 + e^T 1 (dP y / dE y ) for the saturated state 
from the P y vs. E y curve. The calculated dielectric constants 
are also nearly the same for both single- and multi-domain 
states (e yy ~ 150) . Our findings are qualitatively in agree- 
ment with the experimental results of Wada et al4i and Park 
et al4! who have measured the electromechanical response 
of BaTiOs multi-domain single crystals at room tempera- 
ture. For a quantitative comparison, a full 3D calculation is 
required. Nevertheless, our 2D calculation captures the es- 
sential physics of these experiments. 

Next, we study the case with T = 273K. At this temper- 
ature, the free energy in Fi has four minima (±0.21, ±0.21) 
Cm~ 2 and (±0.21, =p0.21) Cm~ 2 , each state corresponding 
to a rhombic phase. As in the earlier case, growth of do- 
mains from the paraelectric phase results in a multi-domain 
state with all four variants, as shown in Fig. 3(a). However, 
for this case, the domain walls are oriented along [10] or [01] 
as dictated by the symmetry of the rhombic phase encoded 
in the anisotropic long-range interaction F es . There is no net 
charge, even for this case. Here also, we apply the electric 
field along the [01] direction. This situation is interesting as 
the [01] direction is not one of the polarization directions at 
this temperature. Because of its superior piezoelectric proper- 
ties, this situation has been studied in numerous experimental 
systems^ii We vary the electric field E y at the same rate as 
in the T = 298A' case. Unlike the earlier case, application 
of the field does not result in transformation to a single do- 
main state. Instead, a stable multi-domain state that has only 
two of the variants (P y > 0) is formed at E y = 5.54 kV/cm 
(Fig. 3b). This configuration is analogous to the engineered 
domains observed in recent experiments^!!. This engineered 
configuration persists till E y ~ 51kV/cm, although individ- 
ual polarization vectors gradually rotate towards the [01] di- 
rection. Figure 3(c) shows the domain pattern at E y = 51.74 
kV/cm. Here we can clearly see that polarization vectors 
have rotated and at the domain wall, the polarization is almost 
aligned with the [01] direction. Thus, the domain boundaries 
serve as a nucleation source for an electric field induced struc- 
tural transition from a rhombic to a rectangular ferroelectric 
phase. The transition is complete at E y = 53.59 kV/cm, as 
is clear from Fig. 3(d) where all the polarization vectors are 
aligned along the [01] direction. If we remove the field, the 
situation shown in Fig. 3(d) remains in a single domain rect- 



angular state. However, if we remove the field before the field 
induced transition, the engineered configuration remains sta- 
ble. 

In Fig. 4(a), we show the variation of (P x ) and (P y ) for the 
situation depicted in Fig. 3. The appropriate single domain 
behavior starting from a polarized state P x = 0.21 Cm -2 
and P y = 0.21 Cm -2 is also shown. Polarization rotation 
and field induced transition at E y ~ 53 kV/cm for the multi- 
domain state is apparent from this figure. Interestingly, the 
transition occurs at a higher field value E y ~ 69 kV/cm for 
the single domain state. This is due to the fact that in the 
single domain state, there are no nucleation mechanisms. Nu- 
cleation sources exist in the multi-domain state due to the do- 
main walls. In the single domain calculations of Bell&, the 
field induced transitions occured at much higher values than 
the experimental values. The present calculation suggests 
that domain walls can help to reduce the field level, similar 
to the effect of dipolar defects in reducing the coercive field 
during switching 10 . Figure 4(b) shows the electromechanical 
response for this situation for both single and multi-domain 
states. These curves are similar to the strain vs. electric 
field curves for experiments where the field is applied along 
non-polar directions^!!!. In the domain engineered regime 
(before the field induced transition), high values of due 
to polarization rotations are found. For example, the multi- 
domain ^33 ~ 567 pC/N and the corresonding single domain 
d 33 ~ 367 pC/N at E y = 49.89 kV/cm are achieved. The cor- 
responding dielectric constants are e yy ~ 750 (multi-domain) 
and e yy ~ 525 (single domain). The difference between the 
single domain and multi-domain electromechanical response 
shows the property enhancement due to the presence of do- 
main walls. 

To conclude, we have used a Ginzburg-Landau formal- 
ism to demonstrate the effect of the underlying domain mi- 
crostructure on electromechanical properties of ferroelectric s. 
To account for nonlocal elastic and electrostatic effects, two 
additional parameters that measure the strength of the long- 
range interactions have been introduced. Our calculations 
show that these long-range parameters are essential to de- 
scribe multi-domain states. In the present work, we have made 
simple choices for these parameters. In principle, these pa- 
rameters should be measured for a given experimental multi- 
domain state for a complete characterization of the material. 
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FIG. 1: Evolution of ferroelectric domains for a simulated (~ 0.85 
/im X0.85 fi m) system at T = 298A". The electric field values 
are E x = 0, E y = (snapshot (a)); E x - 0,E y = 5.54 kV/cm 
(snapshot (b)); E x = 0, E y = 7.39 kV/cm (snapshot (c)); E x = 
0, £ H = 9.24 kV/cm (snapshot (d)); 
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FIG. 2: (a) Variation of (P x ) (red line with circles) and (P y ) (blue 
line with crosses) with the applied field E y for the evolution shown 
in Fig. 1. Also shown is the variation of P x (red solid line) and 
P y (blue solid line) obtained by minimizing only the free energy Fi 
(A = jj, — 0), corresponding to a single domain state, (b) Variation 
of the uniaxial strain (rjyy{E y )) — (f]yy(^)) ( re d line with circles) 
with the electric field E y . The corresponding single domain response 
is also shown (solid red line). 
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FIG. 3: Analogous to Fig. 1 but for T = 273K. The electric field 
values are E x = 0, E y = (snapshot (a)); E x = 0, E y = 5.54 
kV/cm (snapshot (b)); E x = 0,E y = 51.74 kV/cm (snapshot (c)); 
E x = 0, E y = 53.59 kV/cm (snapshot (d)); 
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FIG. 4: (a) Variation of (P x ) (red line with circles) and (P y ) (blue 
line with crosses) with the applied field E y for the evolution shown 
in Fig. 3. Also shown is the variation of P x (red solid line) and P y 
(blue solid line) for a single domain, (b) Variation of the uniaxial 
strain {rjy y (E y )) — (r]yy(0)} (red line with circles) with the electric 
field E y . The corresponding single domain response is also shown 
(solid red line). 



